1 Local MDS for Zero Digits

As in Assignment 2, we will consider only the zeros in the ZIP data set.

# Utility function to plot images. 
plot.zip <- function(x,use.first=FALSE, ...){
  x<-as.numeric(x)
  if (use.first){
    x.mat <- matrix(x,16,16)
  }else{
    x.mat <- matrix(x[-1],16,16)
  }
  image(1:16,1:16,x.mat[,16:1],
        col=gray(seq(1,0,l=12)))
  if (!use.first) title(x[1])
}

zip.train <- read.table("../../unsup/zip.train")
zip.zeros <- zip.train %>% dplyr::filter(V1 == 0) %>% dplyr::select(-1)
dim(zip.zeros)
#> [1] 1194  256
delta <- dist(zip.zeros)

1.1 2-dimensional Configuration of the Data

Using parameters \(k = 5\) and \(\tau = 0.05\) in the function stops::lmds, a 2-dimensional configuration of the data is found.

q <- 2
k <- 5
tau <- 0.05
lmbds.res <- lmds(as.matrix(delta), k = k, tau = tau, ndim = q, itmax = 700, verbose = 10)
#> [1] "niter=100 stress=-251754.80998"
#> [1] "niter=200 stress=-265539.92165"
#> [1] "niter=300 stress=-269891.73297"
#> [1] "niter=400 stress=-271730.14941"
#> [1] "niter=500 stress=-272742.21209"
#> [1] "niter=600 stress=-273206.77601"
#> [1] "niter=700 stress=-273634.96777"
conf <- lmbds.res$conf

The scatterplot is shown below. 9 points are selected in order to try to “cover” the variability in the scatter plot. The points we have chosen are shown in red. They are chosen from left to right along the first dimension and from bottom to top along the second dimension.

# Code to select points. 
find.nearest <- function(new.pos, data){
  # Compute distance to points and select nearest index.
  return(which.min(colSums((t(data) - new.pos)^2)))
}


chosen.indices.localMDS <- c(find.nearest(c(-100,0), conf), 
                                find.nearest(c(-40,0), conf), 
                                find.nearest(c(0,0), conf), 
                                find.nearest(c(20,0), conf), 
                                find.nearest(c(60,0), conf), 
                                find.nearest(c(40,-25), conf), 
                                find.nearest(c(40,0), conf), 
                                find.nearest(c(40,20), conf),
                                find.nearest(c(40,50), conf))


chosen.points.localMDS <- data.matrix(conf[chosen.indices.localMDS,])

plot(conf, as = 1, main=paste0("Local MDS, k=",k,", tau=",tau), xlab = "Dim1", ylab = "Dim2")
points(chosen.points.localMDS, col = 2, pch = 16)

The ZERO digits corresponding to each of the points are plotted below.

par(mfrow = c(3,3))
apply(zip.zeros[chosen.indices.localMDS,],1,plot.zip, use.first = T)

#> NULL

The first five zeros (row-wise) travel from left to right in the first dimension of the scatter plot gained from local MDS. It seems like the first coordinate describes zeros from narrow to wide (in the west-east direction of the drawings). The last four zeros that are plotted show the movement from small values to larger values in the second dimension. It seems like the second coordinates describes zeros from thin thick (in the south-north direction of the drawings). Thus, the first dimensions describe the changing sizes of the zeros, while the second dimensions describes the different line thicknesses of the zeros. This can be seen more clearly when plotting more points, but we have only shown 9 here.

zip.PC <- princomp(zip.zeros[,-1])
scores <- zip.PC$scores[,1:3]

smooth.pc1 <- gam(scores[,1] ~ s(conf[,1],conf[,2]))
smooth.pc2 <- gam(scores[,2] ~ s(conf[,1],conf[,2]))
smooth.pc3 <- gam(scores[,3] ~ s(conf[,1],conf[,2]))

s.pc1 <- smooth.pc1$fitted.values
s.pc2 <- smooth.pc2$fitted.values
s.pc3 <- smooth.pc3$fitted.values

points3D(conf[,1], conf[,2], s.pc1)

points3D(conf[,1], conf[,2], s.pc2)

points3D(conf[,1], conf[,2], s.pc3)

points3D(scores[,1], scores[,2], scores[,3], pch=19,cex=.4, col = 1, 
       xlim = range(scores[,1]), ylim = range(scores[,2]), zlim = range(scores[,3]),
       xlab = "PC1", ylab = "PC2", zlab ="PC3")
points3D(s.pc1, s.pc2, s.pc3, col = 2, add = T)

#surf3D(as.matrix(s.pc1), as.matrix(s.pc2), as.matrix(s.pc3), colvar = as.matrix(s.pc3))

2 ISOMAP for Zero Digits

The same tasks are solved with ISOMAP instead, in order to compare to the results obtained when using local MDS.

k.iso <- 5
ismp <- isomap(delta, k = k.iso, ndim = q)
ismp.points <- ismp$points

The two-dimensional configuration found using ISOMAP can be seen below.

aux.plot <- plot(ismp,n.col=3,main="Output of ISOMAP Algorithm")
points(aux.plot,"sites",pch=19,cex=.6)

chosen.indices.ISO <- c(find.nearest(c(-100,0), ismp.points), 
                                find.nearest(c(-40,0), ismp.points), 
                                find.nearest(c(-10,0), ismp.points), 
                                find.nearest(c(10,0), ismp.points), 
                                find.nearest(c(40,0), ismp.points), 
                                find.nearest(c(25,-45), ismp.points), 
                                find.nearest(c(25,-20), ismp.points), 
                                find.nearest(c(25,0), ismp.points),
                                find.nearest(c(25,25), ismp.points))

chosen.points.ISO<- data.matrix(ismp.points[chosen.indices.ISO,])

plot(ismp.points, main=paste0("ISOMAP, k=",k.iso), xlab = "Dim1", ylab = "Dim2")
points(chosen.points.ISO, col = 2, pch = 16)

After selecting points the red points shown above, the corresponding zeros are plotted in a similar fashion as for local MDS.

par(mfrow = c(3,3))
apply(zip.zeros[chosen.indices.ISO, ],1,plot.zip, use.first = T)

#> NULL
par(mfrow = c(1,1))

As earlier, the first five images show the traversion from west to east in the first dimension, while the last four images show the traversion from south to north in the second dimension. The first dimension seems to describe the … and the second dimension seems to describe the …

TOLKNINGEN HER MANGLER!

3 Select the Tuning Parameters for Zero Digits

3.1 ISOMAP

The local continuity meta criteria is used to select the tuning parameter \(k\) in ISOMAP. The first cell of code is used to calculate the adjusted local continuity meta criteria (function taken from the lecture material).

LCMC <- function(D1,D2,Kp){
  D1 <- as.matrix(D1)
  D2 <- as.matrix(D2)
  n <- dim(D1)[1]
  N.Kp.i <- numeric(n)
  for (i in 1:n){
    N1.i <- sort.int(D1[i,],index.return = TRUE)$ix[1:Kp]
    N2.i <- sort.int(D2[i,],index.return = TRUE)$ix[1:Kp]
    N.Kp.i[i] <- length(intersect(N1.i, N2.i))
  }
  N.Kp<-mean(N.Kp.i)
  M.Kp.adj <- N.Kp/Kp - Kp/(n-1)
  
  return(list(N.Kp.i=N.Kp.i, M.Kp.adj=M.Kp.adj))
}

Next, the function is used to find optimal \(k\). We fix the dimension \(q = 2\). The search for the optimal \(k\) is restricted to integers between 3 and 10. Note that we begin at \(k = 3\) because a lower \(k\) gives fragmented data in the ISOMAP.

Kp <- 5
k.search <- 3:10
nk <- length(k.search)
LC <- rep(NA, nk)
ISOMAP.k <- vector("list",nk)

for (i in 1:nk){
  ISOMAP.k[[i]] <- isomap(delta, ndim=q, 
                            k = k.search[i])
  D2.k <- dist(ISOMAP.k[[i]]$points[,1:q])
  LC[i] <- LCMC(delta,D2.k,Kp)$M.Kp.adj
}

i.max <- which.max(LC)
k.max <- k.search[i.max]
ISOMAP.max <- ISOMAP.k[[i.max]]

plot(k.search, LC, type="b", main=paste0("k.max=",round(k.max,4)))
abline(v=k.max,col=2)

The plot above shows that the optimal \(k\) is \(k =\) 3. The graphical description of the low dimensional configuration corresponding to the optimal \(k\) is given in the following.

ismp.opt <- isomap(delta, k = k.max, ndim = q)
ismp.opt.points <- ismp.opt$points

aux.plot.opt <- plot(ismp.opt,n.col=3,main=paste0("Output of ISOMAP Algorithm with Optimal k=",k.max))
points(aux.plot.opt,"sites",pch=19,cex=.6)

A similar description for 9 chosen digits can be done here. HER MANGLER DET NOE

3.2 Local MDS

The local continuity meta criteria can also be used to select the optimal parameters \(\tau\) and \(k\) in Local MDS. We restrict the search of \(k\) to the integers \(5, 10, 15\) and the search for \(\tau\) is restricted to \(0.1, 0.5, 1\).

Kp <- 10
K <- c(5,10,15)
tau <- c(.1,.5,1)

LC <- matrix(0,nrow=length(K),ncol=length(tau))
LocalMDS.k.tau <- array(vector("list",1),dim=dim(LC))

for (i in 1:length(K)){
  for (j in 1:length(tau)){
    LocalMDS.k.tau[[i,j]] <- lmds(as.matrix(delta), k = i, tau = j, ndim = q, itmax = 100)$conf
    D2.k.tau <- dist(LocalMDS.k.tau[[i,j]])
    LC[i,j] <- LCMC(delta,D2.k.tau,Kp)$M.Kp.adj
  }
}

ij.max <- arrayInd(which.max(LC),.dim=dim(LC))
k.max <- K[ij.max[1]] 
tau.max <- tau[ij.max[2]] 
LocalMDS.max <- LocalMDS.k.tau[[ij.max[1],ij.max[2]]]

print(paste0("k.max=",k.max,"; tau.max=",tau.max))
#> [1] "k.max=5; tau.max=0.5"

HER MANGLER DET NOE

4 t-SNE for Zero Digits.

The same tasks are solved with t-SNE, in order to compare to the results obtained earlier.

tsne <- Rtsne(zip.zeros, dims = q, perplexity = 40, theta = 0)
tsne.points <- tsne$Y
plot(tsne.points, main="t-SNE", xlab = "Dim1", ylab = "Dim2")

plot(tsne.points, as = 1, main="t-SNE", xlab = "Dim1", ylab = "Dim2")

We choose some points in order to “cover” the variance in the two-dimensional representation computed by t-SNE.

chosen.indices.tsne <- c(find.nearest(c(-60,0), tsne.points), 
                                find.nearest(c(-30,0), tsne.points), 
                                find.nearest(c(-10,0), tsne.points), 
                                find.nearest(c(20,0), tsne.points), 
                                find.nearest(c(60,0), tsne.points), 
                                find.nearest(c(-30,-30), tsne.points), 
                                find.nearest(c(-30,0), tsne.points), 
                                find.nearest(c(-30,15), tsne.points),
                                find.nearest(c(-30,40), tsne.points))

chosen.points.tsne <- data.matrix(tsne.points[chosen.indices.tsne,])

plot(tsne.points, main="t-SNE", xlab = "Dim1", ylab = "Dim2")
points(chosen.points.tsne, col = 2, pch = 16)

The corresponding zeros to the selected points are plotted below.

par(mfrow = c(3,3))
apply(zip.zeros[chosen.indices.tsne, ],1,plot.zip, use.first = T)

#> NULL
par(mfrow = c(1,3))

As earlier, the first five images show the traversion from west to east in the first dimension, while the last four images show the traversion from south to north in the second dimension. The first dimension seems to describe the … and the second dimension seems to describe the …

5 Select Tuning Parameter for Zero Digits

The local continuity meta criteria is used to select the tuning parameter perplexity in t-SNE. We are still using the fixed lower dimension \(q = 2\) and we are setting theta \(= 0\). The search for the optimal value of perplexity is constrained to the values c(20, 40, 60).

perplex <- c(20,40,60)
Kp <- 10

LC <- numeric(length(perplex))
Rtsne.k <- vector("list",length(perplex))

for (i in 1:length(perplex)){
    Rtsne.k[[i]] <- Rtsne(zip.zeros, perplexity=perplex[i], dims=q,
                          theta=0, pca=FALSE, max_iter = 500)
    D2.k <- dist(Rtsne.k[[i]]$Y)
    LC[i] <- LCMC(delta,D2.k,Kp)$M.Kp.adj
}

i.max <- which.max(LC)
perplexity.max <- perplex[i.max[1]] 
Rtsne.max <- Rtsne.k[[i.max]]

plot(perplex,LC, main=paste0("perplexity.max=",perplexity.max))
abline(v=perplex[i.max],col=2)

HER MANGLER DET NOE